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A powerful existing technique for evaluating statistical mechanical quantities in two-dimensional 
Ising models is based on constructing a matrix representing the nearest neighbor spin couplings and 
then evaluating the Pfaffian of the matrix. Utilizing this technique and other more recent devel- 
opments in evaluating elements of inverse matrices and exact sampling, a method and computer 
code for studying two-dimensional Ising models is developed. The formulation of this method is 
convenient and fast for computing the partition function and spin correlations. It is also useful for 
exact sampling, where configurations are directly generated with probability given by the Boltz- 
mann distribution. These methods apply to Ising model samples with arbitrary nearest-neighbor 
couplings and can also be applied to general dimer models. Example results of computations are 
described, including comparisons with analytic results for the ferromagnetic Ising model, and timing 
information is provided. 



I. INTRODUCTION 

Just over 50 years ago, Kasteleyn [T] and Fisher and 
Temperley [5] presented analytic combinatorial methods 
for counting dimer packings on a lattice; these techniques 
were soon applied [31 0] to computing the partition func- 
tion for the pure ferromagnetic two-dimensional Ising 
model |6]. These methods continue to be extended 
and improved to study two-dimensional models in statis- 
tical mechanics. Such methods were then extended [71 
to numerically compute the thermodynamic properties 
of disordered magnets. They have allowed for a precise 
and extensive study of the statistical mechanics of dis- 
ordered models (BHTO]. This paper presents a detailed 
description of a numerical approach to implement these 
combinatorial techniques. 

To review the power of these techniques in more detail, 
consider an Ising model where spins on a planar lattice 
can take on one of two values and the energy is given by 
the sum over possible ferromagnetic or antiferromagnetic 
interactions between pairs of neighboring spins. Directly 
evaluating the partition function of this model with N 
spins involves a sum of Boltzmann factors over the 2 N 
spin configurations. But combinatorial techniques allow 
for a much more compact evaluation. The Ising con- 
figurations can be put into correspondence with dimer 
coverings on a related lattice, where dimer coverings are 
choices of edges so that each node of the lattice is in ex- 
actly one chosen edge. Weighted sums Zd over all dimer 
coverings give the partition function Z for the Ising prob- 
lem, with Z = Zd- The sum over all dimer coverings can 
in turn be expressed as the Pfaffian [S] of a weighted 
and signed adjacency-like matrix, the Kasteleyn matrix 
PQ; for a skew-symmetric matrix, the square of its Pfaf- 
fian is equal to its determinant. For the specific case of 
regular lattices and interactions, the Pfaffians can even 
be evaluated analytically by direct diagonalization of the 
Kasteleyn matrix [3j [4] , allowing for exact evaluation of 
thermodynamic quantities and studies of phase transi- 
tions. Pfaffians (or also determinants) ofmxm matrices 



can be defined directly as the sum over permutations 
whose number grows exponentially with m, but the ma- 
trix can be be simplified by column and row eliminations, 
so that the Pfaffian (or determinant) can be evaluated in 
time polynomial in to. As the Kasteleyn matrix used 
here is of size AN x AN, Z can be found for general pla- 
nar Ising models in time polynomial in N . (Note that as 
the arithmetic precision needed for a stable calculation 
of Z depends on ft, there is a ft dependent prefactor for 
the running time which scales roughly as ft 

This technique was subsequently generalized to include 
inhomogeneous couplings between nearest neighboring 
Ising spins by Saul and Kardar [7] for numerical work. 
In the Ising spin glass, the nearest neighbor interactions 
can be either ferromagnetic or antiferromagnetic. For a 
given random choice of couplings of arbitrary sign, the 
Pfaffian of the Kasteleyn matrix can be computed and 
used to derive thermodynamic potentials and suscepti- 
bilities. If the couplings are of fixed magnitude J but 
random sign (the bimodal distribution), the exact de- 
pendence of Z on inverse temperature ft can be written 
as a polynomial in e~ 2 ^ J [7]. More generally, this nu- 
merical approach has allowed for a detailed study of the 
thermodynamics of the 2D Ising spin glass, even with 
continuous disorder distributions. One example of the 
more important recent developments in these algorithms 
has been the application of nested dissection and inte- 
ger arithmetic [5j for computing Z(T) for the bimodal 
distribution in larger systems. Applying Wilson's dimer 
sampling technique j!2j , this numerically exact approach 
has also been used to generate random samples of con- 
figurations of random Ising models [TT]. This sampling 
method bypasses the long equilibration times that arise 
in Markov Chain Monte Carlo methods. 

In this paper, we describe a version of nested dissection 
as applied to the Pfaffian techniques, with the goal of sim- 
plifying and extending the calculation of correlation func- 
tions and sample configurations. For a two-dimensional 
Ising sample with arbitrary nearest-neighbor couplings, 
we describe how this technique can be used to compute 
the partition function, to calculate correlation functions, 
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and to randomly choose sample spin configurations. We 
find that, due to near cancellations of intermediate sums, 
multiple precision floating point arithmetic is needed to 
find accurate results at temperatures of interest for sys- 
tem sizes of size about 20 2 or larger (though the sets 
of integer fields used in Ref. |8] could also be used for 
the partition function in the bimodal case). We note 
that correlation functions were computed at T — by 
Blackman and Poulter for the bimodal case [13] using a 
different approach. We simplify the sampling technique 
used in our previous work by simplifying the matri- 
ces and by using a different approach to maintain com- 
puted correlation functions as spins are sampled. Many 
of these improvements are based on the FIND (fast in- 
verse using nested dissection) algorithm [14] which com- 
putes desired elements of a matrix inverse quickly and 
was developed to compute nonequilibrium Green's func- 
tion applications in nanodevices. This particular flavor 
of hierarchical decomposition is very well suited to the 
geometry of the mapping between two-dimensional Ising 
models and dimer coverings. While much of this algo- 
rithm is implicit in previous work, we assemble these 
methods into a form adapted to studying the statistical 
mechanics of the Ising model, with novel applications to 
computing correlation functions, and emphasize the na- 
ture of the algorithms as a renormalization procedure and 
clarify the sampling procedure. This formulation is also 
significantly faster in practice. We present comparisons 
with analytic results, sample results for the spin glass 
case, and empirical results for the timings. A version 
of the computer code for computing partition functions, 
written in CH — h, is available in the supplemental materi- 
als for this paper at [publisher URL] or by download [15] . 
Extensions of this version of the code have been checked 
against analytic predictions for correlation functions in 
the ferromagnet and against other exact codes for small 
spin glass samples. This code can be used to study pure, 
random bond, and spin glass models. 



II. ISING MODEL, DIMERS, & PFAFFIAN 

In this section, we state the standard Ising spin glass 
Hamiltonian and recall the mapping between Ising spin 
configurations and dimer coverings Qj5] . We also 
review the definition of the Pfaffian of the Kasteleyn ma- 
trix and its relation to the partition function of the Ising 
model. 

A state S of the Ising model in two dimensions on a 
rectangular sample composed of L x x L y spin variables Si 
is given by a choice for each Si , where each Si is restricted 
to Si = ±1. The n = L x x L y sites i lie on a square 
grid. There are 2™ possible spin configurations in the 
state space S. The statistical mechanics of this model is 
governed by the standard Hamiltonian 

(ij) 



where the sample-dependent bond strengths Jij are 
quenched, i.e., fixed in time, and connect nearest neigh- 
bor spin pairs (ij). The spins lie on the nodes of a graph 
G whose edges (ij) connect these nearest neighbor pairs. 
For open or free boundary conditions or for fixed spins 
on the boundaries, these pairs form the edges of a pla- 
nar graph. For periodic boundary conditions, nearest 
neighbor pairs (ij) include bonds that wrap the sample 
around a torus by connecting the top row of the array to 
the bottom row and the right column to the left column. 
In equilibrium at temperature T — , the probability 
P(S) of a spin state S is P(S) = exp-W-lS) / Zj where 
the partition function is Z = X^seS exp - ^^ 5 '. Numer- 
ical derivatives of Z(T) with respect to T allow for the 
computation of energy E(T), the entropy S(T), and the 
heat capacity C(T). Exact sampling will be taken to 
mean that configurations S are generated with the cor- 
rect probability P(S), within numerical accuracy. The 
correlation functions that will be computed by the algo- 
rithm are spin-spin correlation functions, (siSj), where 
the average is taken over all configurations weighted by 
their equilibrium probability, i.e., (siSj) — P(S)siSj. 
Computing these correlations allows for a direct measure 
of the correlation length and the density of relative do- 
main walls. Though we describe the techniques using 
square lattice samples with open boundaries or with pe- 
riodic boundaries, the techniques presented generically 
apply to arbitrary graphs on low-genus surfaces [TJ [5] • 

A given spin configuration S in the Ising model can be 
represented by a set of relative domain walls and by the 
value of a single spin. These domain walls can be defined 
relative to any reference spin configuration S r ; one sim- 
ple choice for S r is the fixed direction configuration S + , 
with all Si = +1. This is the choice that we will use in 
this paper. (Another example choice would be a ground 
state configuration S gs that minimizes %.) The domain 
walls divide the spins into connected sets of spins that 
are either all aligned with or all opposite to the spins in 
S r . These domain walls can be drawn as loops on the 
dual graph Gd- The graph Go has nodes at the cen- 
ter of each (square) plaquette of G. The edges of Gd 
are dual to the edges in G: they are in one-to-one cor- 
respondence, with each edge in Gd crossing one edge in 
G. Given an arbitrary spin configuration S and a nearest 
neighbor pair of spins (ij), the dual edge that crosses the 
bond connecting i to j is in a domain wall if SiSj ^ s£sj. 
For the choice S r = S + , the domain walls separate up 
spins from down spins. As the domain walls are closed 
loops, an even number of domain wall segments meet at 
each node in Gd- 

The configurations in the Ising model may be put into 
correspondence with a complete dimer covering problem 
on a decorated dual graph G* D . For the case we are con- 
sidering, where G is a square grid, each node of Gd can 
be replaced by a Kasteleyn city [3], which is a subgraph 
composed of four fully connected nodes. (Note that a 
Kasteleyn city can be found from a Fisher city [I] by 
Pfaffian elimination.) By replacing each lattice point in 



the dual Gd with a Kasteleyn city, one arrives at the 
decorated dual graph G* D shown in Fig. [I] This larger 
graph allows for a correspondence between domain walls, 
equivalent to Ising spin configurations up to a global spin 
flip Si — > —Si, and dimer matchings on G* D . For each set 
of domain walls, there is at least one corresponding dimer 
covering on the decorated dual lattice G* D . 

The computation of the partition function Z for the 
Ising model, a sum over all assignments of Ising spins, 
can be directly expressed as a related sum over complete 
coverings of either G* , the decoration of the graph G by 
Kasteleyn cities [3] , or coverings of G* D . For sampling and 
computing correlation functions, though, it is simpler to 
start with the decorated dual graph [IIJ [17] . In a pure 
Ising model, summing over matchings on G* corresponds 
to a high temperature expansion |18j . while sums over 
G* D correspond to a low temperature expansion. There 
is a simple correspondence between domain wall loops in 
G* D and spin configurations: given a set of domain walls, 
spins are found by setting spins within a single connected 
region to the same value. 

The configurations contributing to the partition func- 
tion sum correspond to terms in the expansion of the 
Pfaffian of the Kasteleyn matrix [T] : to describe this cor- 
respondence, we first need to define the Kasteleyn ma- 
trix and the Pfaffian sum. The Kasteleyn matrix K is 
a skew-symmetric matrix with non-zero entries for each 
edge of the decorated graph G* D \ the rows and columns of 
the matrix are indexed by the vertices of the decorated 
graph. Skew symmetry implies K a ^ = —K^ a . A non- 
zero entry K a b corresponds to an edge in G* D connecting 
vertices a and b. It is a matrix of size AN x AN. The 
values of the matrix are ±1 for edges internal to Kaste- 
leyn cities. Edges that connect cities have weights with 
absolute value \K a i,\ — exp(— 2/3J^), where spins i and j 
have coupling and the edge ij in G crosses the edge 
ab in G* D . The sign of each K a b is determined by a Pfaf- 
fian orientation of the dimer graph (see, e.g., Ref.[5]). For 
the graph G* D , a simple Pfaffian orientation is that hor- 
izontal edges between Kasteleyn cities are oriented from 
left to right and vertical intercity edges are oriented from 
bottom to top, so that if a is to the left of 6 or 6 is above 
a, K a b > 0. The orientation of edges internal to a Kaste- 
leyn city can then be set as in Ref. [3] or as described in 
Sec. [3] Given a proper choice of signs for K a ^, the parti- 
tition function for the original spin problem on a planar 
graph (without periodic boundaries) can then be shown 
[HE] to be equal to Pf{K), 

Z = Y / e~ mS) =Pi(K), (2) 

s 

where the Pfaffian of K is defined by a sum over permu- 
tations P of node indices, 

Pf(K) =Y,<P)K klh K k2h ---K kmlm , (3) 
p 

with e(P) giving the sign of the permutation P = 
(kx,h, k m , l m ) of the M indices for the nodes of G* D 
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Figure 1: (color online). Correspondence between spin state 
configurations and complete dimer coverings on the decorated 
dual graph G* D for a periodic spin lattice of size L x x L y = 
6x3. A sample configuration of Ising spins Si = ±1 are repre- 
sented by the arrows inside the large circles. These spins are 
coupled by horizontal and vertical bonds of strength Jij . The 
(red) dashed lines indicate the bonds for one example spin. 
The spins and bonds are the vertices and edges, respectively, 
of the Ising model graph G. The nodes of the decorated dual 
graph G* D are drawn as small circles and the edges are in- 
dicated by the thin and thick solid lines. The edges of G* D 
are either internal to a Kasteleyn city (the sets of 4 fully con- 
nected nodes) or connect neighboring cities. Those that are 
internal to a city have a weight of f while those connecting 
cities have a weight w = exp(— 2/3 Jij), where the Jij is the 
coupling strength of the bond crossing the dual edge. An ex- 
ample of a complete dimer covering M corresponding to the 
displayed spin configuration is indicated by the heavy lines: 
such a choice of edges includes all nodes in G* D exactly once. 
The intercity edges belonging to the covering M separate the 
up spins from the down spins and so compose the relative 
domain walls (here we are assuming that the reference con- 
figuration S r is the configuration with all spins up, si = +1). 
Note that for any Kasteleyn city surrounded by 4 spins of 
identical sign, there are 3 ways to arrange the dimers on that 
city. Two examples of these arrangements can be seen in the 
lower left and lower right cities. In other cases, the choice of 
Kasteleyn city edges is uniquely determined by the domain 
walls. 



with m = M/2 = 2L x L y , and the sum is restricted to the 
permutations satisfying the orderings k\ < &2 < • • • < k m 
and k\ < l\, ki < I2, k m < l m . This choice of signs 
forces all domain walls relative to a reference configura- 
tion to enter with a positive sign; it also leads to the 
cancellation of terms such that the many-to-one corre- 
spondence between dimer coverings and spin configura- 
tions becomes one-to-one [TI5] ■ The permutation P of 
indices that enters into the sum represents "matchings" 
or "dimer coverings", i.e., choices of edges e, e\ = {k\,l\), 
. . ., e m = (k m , l m ), such that each node in the dual dec- 
orated lattice belongs to exactly one edge. For proofs of 
the correctness of this mapping see, for example, Kaste- 
leyn's papers [3] and textbook treatments [B]. Note 
that the partition function for a graph of high genus (e.g., 
the three-dimensional Ising model) is impractical to com- 
pute, as it requires a sum over a number of Pfaffians that 
is exponential in the genus |TS] . 

The computation of the partition function on planar 
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graphs is simply given by the evaluation of a single Pfaf- 
fian. Computations on a periodic graph are more compli- 
cated. Kasteleyn described Q] how to compute the par- 
tition functions for dimers on a periodic, i.e., toroidal, 
lattice. Four Pfafhans are computed for four variations 
of K, namely K ++ , K~ + , A" H , and K . We will 
refer to the set of these four matrices by the notation 
A ±=b = {K ++ ,K-+,K+-,K—}. For a given choice for 
r G {+, — } and s G {+, — }, the matrix K rs has matrix 
elements K ab defined according to the standard Pfaffian 
orientation, except for those elements which have end- 
points (a, 6) at opposite ends of the square array: these 
elements correspond to the intercity edges that wrap 
around the graph, leading to a periodic topology. In the 
matrix K rs , if an edge connects the a node a in a city 
that is in column L x — 1 to a node b in column 0, its 
sign is given by r, while if an edge connects a node in a 
city in row L y — 1 to a node for a city in row 0, it has 
sign s. These matrices ii" ±± can be used to compute the 
partition functions Z af> for a = P,AP and /3 = P,AP, 
where P indicates periodic boundary conditions along an 
axis and AP indicates antiperiodic boundary conditions 
[negation of the Jij for all horizontal (vertical) edges in 
a vertical (horizontal) line]. In particular, the partition 
functions are given by linear combinations 



Me(±±) 



(4) 



for a 4 x 4 matrix L [TT [TTj. 
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(5) 



Though the Pfaffian is formally written as a sum over 
a number of permutations that has a number of terms 
roughly exponential in N, the Pfaffian of a general NxN 
matrix can be evaluated in time polynomial in the num- 
ber of nodes, in a fashion similar to computing the de- 
terminant. However, given the two-dimensional nature 
of the graph underlying the matrix A, Pfaffians (or de- 
terminants) and correlation functions can be computed 
much more quickly (a lower power of N) than for a gen- 
eral matrix by splitting the set of nodes geometrically in 
a hierarchical manner |20j . 



III. CLUSTER MATRICES AND THEIR 
OPERATIONS 

In this section, we first give an introductory outline to 
the numerical methods we have implemented for rapidly 
evaluating the partition function and correlation func- 
tions. The details are then described in the subsections 
Sec. |III A| through Sec. |III E| The algorithm for sampling 
configurations is described in Sec. IV 



The introduction to these methods requires the defini- 
tion of the intermediate mathematical objects used, the 
core mathematical steps applied to these objects, and the 
overall organization of these steps to find the Pfaffian for 
the whole sample (or ratios of Pfaffians for correlation 
functions). 

The Pfaffian for the whole sample is computed by com- 
bining information from smaller regions. We can select a 
region A on the decorated dual lattice G* D by choosing a 
loop of Ising spins on the spin lattice G: those nodes in 
G* D that are "inside" the loop (generally the smaller set 
of nodes) will compose the interior set A while those out- 
side the loop compose the exterior, complementary, set 
A. These geometrical regions or clusters have associated 
matrices and factors. The central mathematical objects 
used in this procedure are antisymmetric "cluster" ma- 
trices U A (J) and U-g{J) [HI [20] which depend both on 
the set of spin couplings J — {Jij} and the region A. 
The dependence of U on the spin couplings J that de- 
fine the given realization of a sample will be implicit in 
the remainder of this paper and so we will write Ua for 
Ua{J)- A given cluster matrix is indexed by the nodes 
of the decorated dual lattice that are on the boundary of 
the clusters: if there are m boundary vertices in the clus- 
ter, the matrix U has dimensions m x m. The boundary 
correlations of dimers (and hence spins on the original 
graph) are directly related to the cluster matrices U by 
a matrix inverse. Also associated with each region A is 
a factor z(A), the "partial Pfaffian". This factor repre- 
sents a multiplicative contribution to the overall partition 
function. It represents a sum over dimer configurations 
on the interior of A. 

The core mathematical steps applied to the cluster ma- 
trices are the collection of cluster matrices for neighbor- 
ing regions into a larger matrix and subsequent elimi- 
nation (contraction) steps applied to this joint matrix. 
These elimination steps remove rows and columns from 
the joint matrix that correspond to nodes that are on the 
boundary of the original neighboring regions but are not 
boundary nodes for the union of the two regions. The 
remaining matrix is then indexed by the boundary nodes 
of the larger, unified region. This removal of nodes is car- 
ried out by Pfaffian elimination, a procedure described in 
Sec. lIII Al and one that is similar to Gaussian elimination. 
This directly implements a sum over the the dimer cov- 
erings over edges that are shared by the adjacent clusters 
and incorporates that sum into the partial Pfaffian fac- 
tor. To collect neighboring regions A, with tjia boundary 
nodes, and B, with nig boundary nodes, a square matrix 
of size (mA + mB) X (wa+hb) is filled with the elements 
of Ua and Ub , in block diagonal form, 



M (W,A,B) 



U A W ab 
-Wj b U B 



(6) 



where the matrix W ab is indexed by the boundary nodes 
of A and B and has nonzero elements when a and b 
are the ends of an intercity edge e ab connecting A to 
B. Partial Pfaffian elimination then removes from the 
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matrix rows and columns that correspond to the nodes 
belonging to separating edges in W, while maintaining 
the overall PfafRan. The matrix resulting from elimina- 
tion will be the cluster matrix Uc for the joined regions 
C = A U B, with the matrix again indexed by the re- 
maining boundary nodes. The matrix Uc has dimension 
rn c = rriA + m B — 2\W\. The two steps together, col- 
lection and Pfaffian elimination, will be referred to as a 
"merger" . 

The methods for evaluating the partition function Z = 
Pf (K) are based on relating the Pfaffian of a region of the 
sample to the Pfaffians defined for subregions: by recur- 
sive application of this relationship, the Pfaffian Pf(_ftT) 
of the whole sample can be computed. At the largest 
scale of this recursion, for example, it turns out that we 
can write 



region C. The cluster matrix U^ can be merged with 
Ub- This sums over the configuration sums internal to 
B and the dimer configurations external to both A and 
B, giving a matrix defined on the boundary of A that 
represents the sum over dimer configurations external to 
A, the cluster matrix U-j. At each stage of this recur- 
sion, the cluster matrices Ua and can then be used 
together to find correlations on the boundary of A. It 
turns out that the sums of signed mergers of these two 
matrices gives the spin-spin correlation functions for the 
Ising spins that lie between A and A. 

For reference and to provide a flavor of the methods, 
we present an outline of the steps for computing the par- 
tition function and correlation functions; more detailed 
descriptions of these steps are given in the subsequent 
subsections: 



Pf(K) = a(A*,B*)z(A*)z(B*)IL eEW x e 



(7) 



where the sample is divided geometrically into two re- 
gions, A* and B* , the z(A) and z(B) factors are "par- 
tial Pfaffians" computed recursively. The factors x e 
result from the Pfaffian elimination steps described in 
Sec. Ill A The prefactor a(A*,B*) = ±1 is determined 
by the sign of how the rows and columns of Ua* and Ub* 
are combined. In turn, we can write, for example, 



z(A*) = a(Ai,A 2 )z(Ai)z(A 2 )Il ee w A x e 



(8) 



where the region A* is decomposed into regions A\ and 
A 2 and Wa is the set of edges that connect these two sets 
of nodes. Note that the parity factors a are not strictly 
needed for computing the partition function in planar 
graphs, as all that matters in that case is the magnitude 
of Pf (K) , but they are needed whenever periodic bound- 
ary conditions are used. For numerical stability, pivoting 
operations that permute the rows and columns are used 
and the choice of pivots may be different for the distinct 

The organization of the cluster matrix mergers is di- 
vided into two stages, the up sweep stage and the down 
sweep stage [14] . In each sweep, matrices representing 
neighboring or enclosing regions are merged. The organi- 
zation of these mergers is set by the recursive geometric 
division of the sample. In the up sweep stage, smaller 
cluster matrices Ua and Ub for neighboring clusters A 
and B are merged to create a cluster matrix Uc for the 
union C = A U B of the two clusters. This information 
sums information over smaller scales into information at 
larger scales. The up sweep stage is sufficient to com- 
pute the partition function of a sample. In the down 
sweep stage, correlations (and configuration samplings) 
can be computed. In this stage, the sum of statistical 
weights of all dimer configurations external to a region 
is used to find the sum of statistical weights external to 
smaller regions. If C is the union of clusters A and B, Uq 
gives the matrix encoding the sum of statistical weights 
external to the region C . This matrix is originally found 
by summing over all dimer configurations external to the 



1. From the bond weights Jjy, generate the weights 
Wij = e~ 2 @ J for all neighboring spins in the lattice. 

2. Generate a binary tree T for the geometric subdivi- 
sion of the decorated dual lattice G* D . Each node of 
the tree contains geometric information for a region 
A, the cluster matrices Ua and U^x, and the par- 
tial Pfaffian factors z(A). All non-leaf nodes of the 
tree have pointers to two children representing ma- 
trices for two subregions of approximately the same 
size. The subdivision is terminated at the scale of 
Kasteleyn cities, which are regions that correspond 
to the leaves of T. 

3. Up sweep: starting from the leaves of T, merge sib- 
ling pairs of cluster matrices (Ua,Ub) and factors 
z(A) and z(B) to compute parent matrices Uc and 
partial Pfaffian factors z(C). 

(a) This merging is initiated by collecting the ma- 
trices Ua and Ub along with edge weights for 
the edges W connecting A and B together into 
a joint matrix M (W, Ua,U b ) (see Eq. jfjj). 

(b) Pfaffian elimination then reduces the matrix 
Mq into a set of factors x e and a smaller ma- 
trix Uc indexed by the boundary of AU B. 

(c) Set z(C) = cr(A, B)z(A)z(B)U eeW x e , where 
a(A, B) = ±1 gives the total parity of 
row/column permutations that were used in 
the rearrangements of M in preparation for 
Pfaffian elimination and the parity of permu- 
tations used for pivoting steps during the Pfaf- 
fian elimination. 

(d) These up sweep steps are carried out recur- 
sively, merging clusters up to, but not includ- 
ing, the last pair A* and B* representing the 
initial division of the whole sample. 

4. The two largest clusters for A* and B* are then 
merged according to the choice of boundary condi- 
tions: 
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(a) For open or fixed boundary conditions, sim- 
ply merge the two top-level cluster matrices 
U A" and Ub*- In this case, Pfaffian elimina- 
tion eliminates all rows and columns and the 
partition function Pf(if) is given by Eq. 

(b) For periodic boundaries, merge the Ua* and 
Ub* along one of the rows or columns sepa- 
rating them (there are either two rows or two 
columns separating them for periodic BCs) 
into a matrix Uc* ■ Then connect the matrix 
Uc* with itself along a remaining row to gen- 
erate two matrices Uq» and Uq» , the former 
using positive weights for the wrapping edges, 
the latter using negative weights. Eliminate 
those connecting edges. Then include wrap- 
ping edges along the remaining axis, again us- 
ing negative and positive edge weights for each 
of the Uq* . The resulting eliminations give 
scalars: these overall weights are the Pfafhans 

Pf(iY ±± ). 

(c) Compute Z p < p , Z AP < P , Z P > AP , Z AP < AP from 
linear combinations of Pf(A" ±=t ), as given by 
Eq. @. 

5. Stop here if only the partition function is required. 
Continue to the next steps to compute correlation 
functions. 

6. Down sweep: descend the tree T, computing clus- 
ter matrices for complementary regions and merg- 
ing interior and exterior matrices to find correlation 
functions: 

(a) Use the results of the up sweep to initialize the 
two top level complementary cluster matrices 
via U^r = Ub* and U-g^ = Ua* ■ 

(b) If periodic boundary conditions are used, 
merge Ua* and U-^- and merge Ub* with Ujp- 
using the four different choices for wrapping 
edge weights, i.e., select (r,s) from (±,±). 
Use these mergers to set up four parallel trees 
for further descent. 

(c) For all down sweep steps for a planar Ising 
model or further descending steps in the case 
of periodic boundary conditions in each of the 
four trees: 

i. Given a parent C with known Uq and 
children A and B, merge the parent ma- 
trix XJ-q with Ub to generate matrices 
for regions complementary to A, as in the 
FIND algorithm [H]. 

ii. Also merge Uq with Ua to generate Ug. 

(d) Compute correlation functions between spins 
on the corners of any given region A by signed 
merging of and Ua- (To find correla- 
tion functions for periodic boundary condi- 
tions, compute the correlation function as 



the weighted sum over four trees as given by 



Eq. (13).) 



A. Pfaffian elimination 

Pfaffian elimination simplifies a matrix by setting cho- 
sen elements in a row to zero while maintaining the Pfaf- 
fian of the matrix as an invariant. This elimination pro- 
ceeds by a process similar to Gaussian elimination for 
general matrices, but is applied to skew-symmetric ma- 
trices |21j . In Gaussian elimination, the lower triangular 
elements are set to zero and the determinant is the prod- 
uct of the diagonal elements. In Pfaffian elimination, the 
diagonal elements of a given skew-symmetric U are zero 
and Pfaffian elimination aims to set all elements that 
are more than one step off of the diagonal to zero. The 
Pfaffian of the matrix is the product of the remaining 
elements in even-indexed rows (given that the first row 
has index 0). 

Pfaffian elimination can be defined inductively for a 
skew-symmetric matrix U. Each step simplifies one row 
to a single non-zero element. Suppose that Pfaffian elim- 
ination has been carried out for rows with index less than 
i, where i is even and the rows are indexed starting with 
row and that the element in row i and column i + 1 
is non-zero. Then multiples of row i and column i + 1 
can be added to rows and columns of higher index to 
zero out the remaining elements of row i and column i. 
This addition of rows and columns simplifies the matrix 
while the Pffafian is unchanged, in the same fashion as 
row and column additions in a matrix do not modify its 
determinant. Specifically, for j > i + 1, column i + 1 is 
multiplied by —Ui.j/Ui^+i and added to column j + 1 and 
row i+1 is multiplied by the same prefactor and added to 
row j + 1 [H] . Note that the odd rows do not contribute 
to the Pfaffian when the elimination in the previous even 
row is completed, so that elimination is applied only to 
even rows. We use pivoting of the rows and columns 
that are to be eliminated to improve numerical stability. 
A pivot is an interchange between indices c and d: the 
elements of row c are swapped with the elements of row 
d at the same time columns c and d are swapped. As we 
use it here, Pfaffian elimination is often carried out only 
for some subset of rows. Note that rows/columns that 
are not to be eliminated are not considered for pivot- 
ing. The permutations due to pivoting operations place 
the element with the largest available magnitude in the 
superdiagonal position, before the elimination is carried 
out. Each pivot leads to a change of sign in Pf({7) which 
is accumulated in the prefactor a. 

Mathematically, Pfaffian elimination carried out for all 
rows can be used as a factorization scheme, similar to LU 
factorization via Gaussian elimination [21]. The Pfaffian 
elimination procedure applies linear operations to U so 
that LU L T = F where L is a lower triangular matrix 
and F is zero except for the superdiagonal elements. The 
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inverse of a skew-symmetric matrix U is then 

XJ- 1 = L T F~ 1 L; (9) 

this procedure of elimination and matrix multiplication 
is used to find matrix inverses in the sampling of Ising 
spin configurations (see Sec. IV and Ref. [H]). In the 
mergings of matrices used here, the rows and columns 
corresponding to nodes on the boundary of the joined re- 
gions are kept, while the rows and columns correspond- 
ing to nodes shared by the joined regions are eliminated. 
The eliminated rows and columns have superdiagonal el- 
ements which are multiplied together to give a partial 
Pfaffian while the rows and columns for the new bound- 
ary are carried onto the next stage. Physically, by elimi- 
nating rows and columns corresponding to nodes internal 
to a geometric region, these steps "integrate out" degrees 
of freedom internal to the new cluster. 



B. Geometric dissection 

Computations for sparse matrices that are derived 
from two-dimensional graphs can be very efficiently car- 
ried out using the important technique of nested dissec- 
tion [20] . The row and column indices of the matrix 
correspond to a numbering of the nodes in a graph. The 
idea behind nested dissection is to hierarchically subdi- 
vide the matrix according to row and column indices that 
index nodes for distinct compact regions. When subdi- 
viding a region into two child regions, the separator for 
this subdivision can be taken to be either nodes that lie 
between the two regions or a set of edges that connects 
the two compact subdivisions. By "compact", we mean 
regions of size N whose boundary scales as 0(s/N). The 
result for the parent region is found by separately com- 
puting the results for the two child regions and stitch- 
ing those two results together using the separator. As a 
separator can be found with 0(y/N) nodes for a matrix 
of scale N x N (i.e., of order L for a spatial region of 
size L 2 = N), with the two subproblems of comparable 
size, the computation at each scale is for matrices of size 
O(VN) x O(VN) [2D]. The work at each scale is there- 
fore much less than for the dense case where the problem 
cannot be efficiently subdivided and one needs to con- 
sider matrices of size N x N. The first application of 
nested dissection to efficiently computing spin glass par- 
tition functions is described in Ref. [5] . The use of the 
general concept of nested dissection for sampling dimer 
configurations was proposed in Ref. [12j and carried out 
for Ising spin glasses in Ref. [IT] . 

In the form of nested dissection [2D] used for dimer 
sampling [12], a set of nodes in a graph is selected as 
the separator. This is the form we previously used [TTJ 
for sampling Ising spin configurations. Here, we instead 
use an edge separator with each separating edge having 
one node in each child region. An example approach 
that inspired our method is the FIND (fast inverse us- 
ing nested dissection) technique, which computes some of 



the elements of an inverse matrix, as used in computing 
non-equilibrium Green's functions in a two-dimensional 
quantum device. The asymptotic run-time of computing 
the Pfaffian with either node or edge separators scales 
with N in the same way, i.e., as 0(iV 3 / 2 ) but the FIND 
approach has several advantages for studying the Ising 
model. These advantages include simplifying the struc- 
ture of the code as well as allowing for more direct com- 
putations of the inverse matrix elements and the Pfaffian 
ratios used to sample configurations. 

In the Ising model, the decorated dual graph G* D for 
an L x x L y square sample with periodic boundaries can 
be recursively divided by splitting it either horizontally 
or vertically at each stage into smaller rectangles. Fig. [2j 
gives an example of this dissection. The geometric dissec- 
tion of the system into smaller rectangles is described by 
a binary tree T. Each rectangle is an array of Kasteleyn 
cities. The leaves of the tree consist of I x I arrays, that 
is, individual Kasteleyn cities, so that the corresponding 
cluster matrix Uy for a city Y is a 4 x 4 matrix. At each 
stage of the dissection, the graph is divided along the 
axis with the shortest length and as close to the middle 
of the rectangle as possible. This division splits the city 
set by cutting the edges which join neighboring cities; 
these cut edges comprise the separating set W at each 
stage. It is important to note that the separator W has 
a corresponding set of Ising spins: the spins that lie be- 
tween the two geometric regions and are separated from 
each other by the edges in the set W (see Fig. |4| . The 
top of the tree T has no boundary and so is associated 
with a null matrix at the end of the algorithm. However, 
for efficiency in collecting partial results, the region C* 
corresponding to the whole sample has matrices associ- 
ated with it during intermediate stages of the calculation. 
The two rectangles A* and B* that result from the first 
division of the sample are the first non-empty regions, 
with 4*UB* = C*. 

This tree structure is used to organize the elimination 
steps in the FIND-based technique, which consists of two 
stages: an up sweep which produces the partition func- 
tion of the system by Pfaffian elimination, and a down 
sweep which may be used to find inverse matrix elements, 
bond probabilities, or correlation functions. These stages 
may be understood as a reorganization to move informa- 
tion about dimer correlations on the region boundaries 
from one scale to another. First, in the up sweep stage, 
the cluster matrices, which represent boundary informa- 
tion about couplings that remains after summing over 
internal degrees of freedom, are joined with cluster matri- 
ces in neighboring regions to generate cluster information 
at a larger scale, for the joint region. This is repeated 
until the aggregate thermodynamic properties of the en- 
tire sample are found. Next, in the down sweep stage, 
this information may be propagated back down to give 
correlation results at the smallest scales, and all scales 
in between. This propagation is effected by merging the 
correlation information exterior to a region with the in- 
terior correlation information. 
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Figure 2: (color online) Depictions of the geometric dissection 
tree T for an example decorated dual graph. The original 
Ising system has L x x L y = 3 x 2 spins on a periodic graph; 
the spins are indicated by gray circles. The decorated dual 
graph has 3x2 Kasteleyn cities (the sets of four fully con- 
nected nodes). At each stage, a parent rectangle of Kasteleyn 
cities is divided into two roughly equal sibling rectangles, the 
children. In the algorithms described here, each rectangle A 
has geometric information, two cluster matrices Ua and Ug, 
and a partial Pfaffian factor z(A). (a) The nested dissection 
in real space. The two largest subregions with boundaries, 
A* and B* , indicated. The region C* in the final stage (not 
shown for clarity) is the union of A* and B* , C* = A* U B*; 
at the end of the algorithm, it has no boundary but a ma- 
trix corresponding to this region is used as a working matrix 
when using periodic boundaries, (b) A diagram of the re- 
sulting nested dissection tree T. The leaves of the tree are 
Kasteleyn cities. The root of the tree has no cluster matrix 
associated with it, as the sample has no boundary, though it 
has a Pfaffian factor associated with it, which is used to find 
the partition function for the whole sample. 



C. Up sweep stage 

The matrix operations for the Pfaffian eliminations 
carried out in the up sweep stage can be illustrated by an 
example of the first steps of the algorithm. Describing 
these first steps allows us to display the matrices used 
and their correspondence to the graphs showing dimer 
correlations. Subsequent steps use larger matrices, but 
have the same structure. 

The lowest level steps merge the cluster matrices for 
two neighboring Kasteleyn cities. Let two such cities be 
denoted by A and B. These cities each correspond to 
neighboring nodes on the dual square lattice Gjj and are 
two of the leaves of the tree representing the geometric 
dissection of G* D . The corresponding Kasteleyn matrices, 
Ua and Ub, are the simplest cluster matrices. The rows 
and columns of A and B are each indexed by four nodes, 
so Ua and Ub are each of size 4x4. In general, due 
to skew symmetry, only the upper triangular portion of 
each cluster matrix need be stored in memory. The ele- 
ments above the diagonal in the matrices Ua and Ub are 
displayed in Fig. |3ja) . Each cluster matrix has a weight 
associated with it, a partial Pfafnan that accumulates 
the weights of eliminated rows, which is initialized to be 
unity, z(A) = z{B) = 1. The entries of each matrix have 
weight of magnitude 1, with signs appropriate for Kaste- 
leyn cities (HJ. For the sign conventions and numbering 
scheme show in Fig. [3]ja), all elements of Ua and Ub 
are positive in the upper triangular section. The edges 



joining cities are directed in the positive x and positive 
y directions, so that the matrix elements are non- 

negative for nodes k in cities to the left of or below the 
city containing node I. Here, the separator W consists 
of a single edge e joining city A to city B. The weight 
of this edge is w = exp(— 2/3 J^ ), where is the bond 
weight on the connection between spins i and j that is 
perpendicular to this dual edge e. Note that this edge is 
not connected to the rest of the lattice and so will be be 
eliminated when merging A and B. To carry out this re- 
duction, the elements of Ua and Ub and the edge weight 
are copied into a joint temporary matrix Mq(W, Ua, Ub), 
which is of size 8x8 (28 upper triangular elements). This 
edge to be eliminated is then placed in the first row of 
the matrix by permuting the rows of the joint matrix 
to obtain Mi(W,Ua,Ub)- Whenever two rows are in- 
terchanged, an overall minus sign is introduced into the 
Pfaffian factors. In this simple case, only one Pfafnan 
elimination is applied to Mi({w}, Ua, Ub)- This elimi- 
nates the connections of the ends of the connecting edge 
to the rest of the boundaries of A and B, giving a ma- 
trix with the first superdiagonal element xi^ is non-zero, 
but the rest of the first row eliminated. The remaining 
rows, the third through the last rows, define the new 
cluster matrix Uc(W, Ua, Ub)- This matrix encodes the 
correlations along the outer boundary of C, the region 
composed of the two joined cities. This contracted ma- 
trix is generally not sparse; see Fig. [3j Using Eq. (|8j), 
the partial Pfaffian factor that is stored along with Uc 
is z(C) = —z(A)z(B)x, where the minus sign is included 
because of the row interchange. 

This process of merging adjacent subgraphs of G* D , 
which uses Pfaffian elimination to remove adjacent 
boundary nodes, is repeated at each scale up to the sys- 
tem size L. Generally, neighboring regions A and B are 
joined together by copying their entries into a joint ma- 
trix Mo, adding the weights of connections for the set of 
n edges W that join A to B (W is indicated by jagged 
lines in Fig. [4]), permuting the joint matrix to give Mi, 
and then eliminating the first \W\ rows. The portion of 
the matrix that is indexed by the boundary of C = AUB 
is the larger scale cluster matrix Uc- The product of the 
superdiagonals on the even eliminated rows are used to 
find z{C) = a(P)z(A)z(B)U 2 i ^ l ~ 1) x ili+1 , where a(P) 
is the sign of the permutations carried out in assembling 
and carrying out pivot eliminations during the elimina- 
tion and the a^i+i for even i for the eliminated edges 
are the superdiagonal elements remaining after Pfaffian 
elimination. This process is an exact real-space renor- 
malization process on the space of cluster matrices. At 
each scale, the cluster matrices represent geometric re- 
gions whose interactions are computed using their adja- 
cent boundaries, though each of these clusters has many 
internal degrees of freedom. At the largest length scale, 
the Pfaffian of the remaining O(L) x O(L) matrix is mul- 
tiplied by the products z(A*) and z{B*) resulting from 
all lower level mergers to gives the Pfaffian of the en- 
tire Kasteleyn matrix, i.e., the partition function Pf(if). 



This procedure of Pfaffian elimination and collection of 
superdiagonal elements preserves the overall Pfaffian at 
each stage, since Pfaffian elimination maintains the Pfaf- 
fian as an invariant and the Pfaffian of a matrix with only 
superdiagonal elements in the odd rows is just the prod- 
uct of those superdiagonal elements. The total number 
of operations in a full up sweep is dominated by the last 
merger and is of order 0{N 3 / 2 ). 



D. Down sweep stage: overview 

While the up sweep stage can be used to compute 
a global quantity, e.g., the partition function for given 
boundary conditions, the subsequent down sweep stage 
provides a powerful method for computing spatial infor- 
mation such as spin-spin correlation functions. Correla- 
tion functions at multiple scales can be computed in a 
single down sweep, while multiple down sweeps are used 
to generate sample configurations (see Sec. IV), 

The down sweep stage descends the geometry tree T, 



(a) 



recursively computing new cluster matrices U-j. These 
matrices contain information about sums over dimer con- 
figurations for the exterior A of the geometric regions A. 
These exterior clusters are merged with the interior clus- 
ter matrices that were computed on the up sweep to find 
spin-spin correlation functions. The computed correla- 
tion functions, i.e., the thermal averages (siSj), are for 
pairs of spins i and j that border a geometric cluster 
A: these spins lie between A and A. In our current im- 
plementation of the down sweep stage, we compute all 
pairwise correlations between the four spins that are on 
the corners of each rectangular region. The results of 
this computation include correlations between all pairs 
of neighboring spins, as these are on the corners of the 
region around a single Kasteleyn city (a 1 x 1 region). 



E. Description of the down sweep stage 

The down sweep stage uses as initial data the cluster 
matrices Ua for each node A of the tree found during the 
up sweep. As in the FIND method [14] . this initial set 
of cluster matrices is then used to calculate cluster ma- 
trices U-£ for the complementary (i.e., exterior) regions 
A. These matrices encode the boundary correlations of 
dimer matchings resulting from summing matchings over 
the portion of the sample surrounding a geometric region 
A. This is to be compared with the cluster matrix Ua 
which contains information about dimer correlations be- 
tween its boundary nodes resulting from summing all of 
the dimers within the region A. 

At the highest level, where there are two regions A* 
and B* , U^r — Ub* and U-j^r — Ua*, up to permuta- 
tions of rows and columns due to differing indexing of 
the boundary nodes, as A* is exterior to B* and B* is 
exterior to A*. The complementary matrix U-j for a re- 
gion A at a lower level is computed by merging XJ-q with 
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Figure 3: Depiction of merger of Kasteleyn cities A and B 
connected by an edge of weight w = exp(—2f3Jij), with w 
set to 0.5. (a) The rows and columns of the Ua and Ub are 
indexed by nodes {0, 1, 2, 3} and the upper triangular part of 
these skew-symmetric matrices is shown with the diagonal of 
zeros omitted. For example, the upper left elements shown 
in the matrices here are in row and column 1. The signs 
of the connections correspond to a Pfaffian orientation [3] to 
consistently count spin configurations, (b) The matrix Mo is 
the result of collecting Ua and Ub and the edge weight w = 
0.5 and is indexed by nodes through 7, with indices from 
A for the initial part {0,1,2,3} and indices from B for the 
second half {4, 5, 6, 7}. Permuting rows/columns and 1 and 
then 1 and 7 gives the matrix Mi. These permutations place 
the nodes for the edge to be eliminated in the first two rows 
of Mi. (c) After using Pfaffian elimination to remove rows 
and columns and 1, one is left with a superdiagonal element 
at (0, 1) of 0.5 (no pivoting is possible in this case) giving a 
partial Pfaffian z(C) = 0.5 and the next generation cluster 
matrix Uc, indexed by the remaining 6 nodes numbered as 
shown. 



Ub, where C is the parent region for the siblings A and 
B. The matrices Uq and Ub are placed into a larger 
matrix and the edge weights for those edges whose ends 
are shared by these two boundaries are included. Those 
edges shared by C and B are eliminated by Pfaffian elim- 
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Figure 4: Example of a higher level step of the up sweep stage. 
The Ising spins on the original lattice whose couplings Jij are 
relevant to the calculations through this step are indicated by 
the large grey circles, while the nodes of the decorated dual 
lattice are indicated by the medium-sized and small circles. 
The solid straight and jagged lines indicate edges belonging 
to the decorated dual graph. Two geometric regions on the 
dual lattice, A and B, each containing 3x3 Kasteleyn cities, 
are denoted by the (blue) dashed squares. The dual edges W 
that separate A and B are drawn as jagged lines and join the 
filled in medium-sized nodes. The small nodes are interior to 
the regions A and B while the larger nodes (medium-sized 
open and filled circles) form the borders of A and B. In the 
up sweep stage (see Sec. Ill C I, given the cluster matrices Ua 
and Ub for A and B, the separating edges belonging to the 
separator W are integrated out via Pfaffian elimination, leav- 
ing a cluster matrix for C that includes the entire subgraph 
shown. The cluster matrix Uc is indexed by its border nodes, 
i.e., the medium-sized open circles connected by the dashed 
(blue) rectangle. 



ination and what remains is the cluster matrix U^, in- 
dexed by the nodes adjacent to A. The entire tree is 
descended in this fashion, thus generating complemen- 
tary matrices and correlations between corner spins for 
each region in the geometry tree T. A step of this process 
is diagrammed in Fig. [5] 

As the U-j are computed, the clusters Ua and U-^ can 
be merged via Pfaffian elimination. By comparing the re- 
sults found using different signs for the connecting edge 
weights, the spin-spin correlations on the original lat- 
tice can be computed. To explain this computation of 
correlations, we continue to suppose that domain walls 
are defined using an all spin up reference configuration 
S r = S + , so that neighboring Ising spins of opposite sign 
are separated by a domain wall. Then the Boltzmann 
weig ht e-^W of a given spin configuration S is equal 
to the product c r H ee x(s) w e of all weights w e of edges 
e that make up the domain wall set M(S, S r ) on the 



dual graph with c r 



, which is a sample and 



reference state dependent constant . Consider two spins 
located at i and j in G. In a given spin configuration 
the spins are separated by either an even number or 
and odd number of domain walls in M. The spins have 
equal orientations, Sj = Sj, if and only if a path in G 
between the two spins crosses an even number of domain 
walls. So the correlation function can be found from the 
average parity of domain walls between the two spins i 




Figure 5: (color online) Diagram of a sample merging of clus- 
ter matrices in the down sweep stage for the regions indicated 
in Fig. [4] The gray and black larger circles indicate the loca- 
tions of the Ising spins in the original square grid. (The black 
spins are the corner spins for region C). The matrix U-q de- 
scribes the dimer correlations between the nodes that touch 
the outer (red) dashed line labeled C. This matrix sums over 
correlations external to the spins shown. This cluster matrix 
V-Q is merged with Ub by eliminating the edges shown by the 
jagged solid lines. The region B is indicated by the dashed 
(blue) square on the right of the diagram. The result of the 
merger is the matrix describing correlations among the 
nodes on the boundary of A, which is shown by the labeled 
square (red) dashed line. 



and j. 

Given a choice of couplings Jij with chosen boundary 
conditions and temperature, let the equilibrium fraction 
of spin configurations S with Sj = Sj (sj 7^ Sj) be given 
by P(si — Sj) [respectively, P(si 7^ Sj)]. Let i — > j 
indicate a path of length £ between i and j built up of 
nearest neighbor pairs (i, k\), (k±, h?), ■ ■ ■ , (ke-±,j). For 
I = 1, the path is just the single bond The partition 

function under the constraint that s, 7^ Sj is 



fSH{S) 



(10) 



and can be represented as the restricted sum over match- 
ings M in G* D 



E 



e(P)Y[w(e). (11) 



{M\i—tj crosses odd # edges in M} 



eeM 



In this formula, the sum over matchings is understood to 
be restricted to edge choices that obey the restrictions 
described below Eq. Q and e(Pjy) is the sign of the 
permutations in the listing of the nodes in the matching 
M, so that the many-to-one mapping of dimer coverings 
to spin configurations is effectively turned into a one-to- 
one mapping by cancellation of oppositely signed terms. 
The restriction is to matchings such that the bonds for 
Ising spin pairs in the path i — > j are crossed by the 
dual edges in the matching M an odd number of times. 
Note that this restriction is independent of the exact 
path i — > j and only depends on the endpoints i and 
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j. A similar correspondence (with the sum over an even 
number of crossings) holds for expressing Z Si — s ., with 



z s 



= z. 



To compute these restricted partition functions, we 
compare Pf(JT) with Pf(Ki-+j), where is a mod- 

ified Kasteleyn matrix. All weights for the edges in G* D 
that cross the chosen path i — > j contained in G are 
negated in this modified matrix. That is, the elements 
of Ki^j are the same as those in K except where an 
edge e € G* D is crossed by the path i —> j: in that case, 
the weight w e in K is replaced by — w e in K^j. These 
negations reverse the signs of the weights of matchings 
M with an odd number of edges of that cross i — > j while 
maintaining the sign of matchings with an even number 
of edges of M crossing that path. To efficiently carry out 
the computation of Pf (K^j), the path i — ¥ j is chosen 
to cross edges that connect a cluster A to its complement 
A. That is the path connect spins that lie between A and 
A. An example showing the negated dual edges is shown 
in Fig. § 

These correspondences allow us to write the spin cor- 
relation function for a planar Ising model in the form 

(SiSj) = P(Si = Sj) - P(Si ^ Sj) 

= Z \Z S .— S . Z s 



1 



Pf(if) 



e e ( p ^) n 



{M\c 



e£M 



e n w ^ 

{M|odd eeM 

Pf[M (W i -+ j ,U A ,UA)] 
W[M (W,U A ,Ux)] ' 



(12) 



The modified list of weights W^j is the set of 
weights with negated values for all dual edges crossed 
by the path (any path) from i to j, i.e., i — s- j, 
where again i — > j lies between A and A. Note 
that Pf(K) = Pf[M (W)]z(A)z(A) and Pf = 
Pf [Mo(Wi^jj]z(A)z(A); the cancellation of the common 
factor z(A)z(A) gives the last step in the above equation. 

This representation of the spin correlations in Eq. ( 12 ) 
defines the procedure for their computation. Correla- 
tions between two spins that lie between a region A and 
and its complement are computed by merging the 
two matrices Ua and once using the original weights 
and again using the modified (partially negated) weights. 
The ratio of the two resulting Pfaffians gives the spin-spin 
correlation value. We note that a different approach has 
been used to compute correlation functions at T — [12] , 
where paths between frustrated plaquettes are the basis 
of the representation in the ground state. For the exam- 
ple shown in Fig. [6j the correlation function is calculated 
between the two spins diagonally opposite (top left and 
bottom right) between the outer and inner set of nodes, 



whose correlations are given by Ua and U-r. The choice of 
signs for W could also be modified to compute multispin 
correlations. 




Figure 6: (color online) Diagram of a calculation of a corre- 
lation between spins i and j. This calculation uses informa- 
tion computed during the both the up sweep and down sweep 
stage. The matrices JJ-q and Uc are merged using the edge 
weights on the jagged lines. These edges form the set W. 
Two mergings are calculated: one for all positive weights for 
the edges in W and one where the weights are negated for the 
thicker (green) edges. This gives two Pfaffians for the whole 
sample. The first Pfafhan has positive contributions from 
configurations with either an odd or even number of dimer 
choices (i.e., domain walls) between spins i and j. This is the 
partition function for the whole sample. The second Pfaffian 
is the difference between the partition function constrained 
to have an even number of domain walls between i and j and 
the partition function constrained to have an odd number of 
domain walls between i and j. 



The calculation of correlations is simplest for planar 
graphs (Ising models with open or fixed boundary con- 
ditions). If periodic boundary conditions are to be used, 
four different mergers of the two top level matrices Ua* 
and U b* are computed. These mergers are computed 
for all possible pairings of negative or positive weights 
for bonds that connect the top row to the bottom row 
of cities or the rightmost column to the leftmost col- 
umn of cities, as justified in Sec. IIIC and Ref. [llj . 



The descent of the tree for each is carried out start- 
ing from each of these four choices. There will then be 

. -t- -t- 

four complementary cluster matrices, LL' , for each ge- 
ometrical region A. Each complementary cluster matrix 
will have its own four partial Pfaffian factors z r ' s (A). 
A spin-spin correlation for periodic boundary conditions 
(as given by K ++ ) is then the ratio of two weighted 
sums. The weighted sum in the denominator is the to- 
tal partition function divided by z(A). The sum in the 
numerator is the same linear combination of the Pfaffi- 
ans but with the weights W negated on edges that cross 
the path i — > j (i.e., using the weights Wi-yj). The 
partition function for periodic boundary conditions is 
Z p > p = | J2( r . s ) Pi ( Kr,s )- The iactor Z ( A ) is common 
to all terms Pf (K) = z(A)z(A) Pf [M ] in the weighted 
sums and so can be cancelled out. This gives the re- 
sult that the spin-spin correlation function on a periodic 
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lattice is the ratio 

, . _ E(r,.) e {(±,±)} Pf[M (WWi, U A , U^)]z^(A) 



(r, S )e{(±,±)} 



Pi[M (W,U A> U* 



(13) 



So the correlation function computations, which require 
the evaluation of two Pfaffians on a planar graph, require 
8 Pfaffians on a torus for each pair of spins. The com- 
putation time for the correlation functions for the corner 
spins on all regions in practice requires about 10 times 
the amount of computing time as finding only the parti- 
tion function Z. 



IV. SAMPLING 

Exact sampling methods select independent configura- 
tions according to their probability in the whole sample 
space. We consider here the problem of generating a 
sample configuration of a system with probability pro- 
portional to the Boltzmann weight e~P n . As a contrast 
with direct sampling, consider Markov chain Monte Carlo 
(MCMC) methods. In an MCMC method, a sequence of 
configurations is generated by randomly chosen updates; 
if the update choices obey detailed balance and can reach 
all possible configurations, in the limit of large times this 
sequence will generate sample configurations from the 
Boltzmann distribution [22] ■ The number of Monte Carlo 
updates needed for the approach to fair sampling is often 
unknown and can be very long. However, MCMC meth- 
ods can generate exact sampling if coupling from the past 
[25] can be used to guarantee fair samples (but not nec- 
essarily fast mixing times). However, no known coupling 
methods are practical for Ising spin glass models at low 
temperatures [24]. Markov chain Monte Carlo methods 
are of course of great practical use, but the availability 
of exact sampling in some cases provides for a very use- 
ful comparison and the potential for much more rapid 
calculations for large glassy systems. 

The direct sampling methods we use [llj to generate 
random Ising spin configurations are based on the map- 
ping between dimer and Ising model configurations and 
on dimer sampling methods [12] that use nested dissec- 
tion. By directly selecting a random matching on the 
decorated dual graph G* D with the proper probability, 
we fairly select a set of relative domain walls and hence 
the relative orientations of the spins on the original lat- 
tice. We report here on a modification of the method 
used in Ref. [TT] ; here we use the edge separators W [H] 
described in Sec. Ill rather than a node separator [12] , 



This modification significantly speeds up the sampling 
algorithm for the Ising model, as the dimension of the 
matrices to be factorized are reduced by a factor of three 
from those used in Ref. [IT] The implementation of the 
algorithm is also simplified. 

As in the computation of the partition function and 
correlation functions, the direct sampling calculations 
rely on a geometric dissection. The tree used for sampling 



differs some from that described in Sec. Ill for computing 
partition functions and correlation functions. For sam- 
pling configurations with periodic boundary conditions, 
we start this modified dissection with a cluster C* which 
is formed by joining the two system halves A* and B* 
along a single line of spins. The cluster C* includes all of 
the nodes in the decorated dual graph G* D , but does not 
include the edges at the top or right that connect the top 
row of nodes to the bottom row or the right column to 
the left column. These edges that are left out are those 
used to complete the periodic boundary conditions. See 
Fig. ga) for a dr awing of C* and the initial separator. 
All of the edges internal to the region C* are contracted 
out by Pfaffian elimination in an up sweep to give the 
cluster matrix Uc- This matrix M(C*) is used in the 
first stage of spin assignment. In this first stage, the Ising 
spins that form the bottom row of the sample are chosen. 
As the probability distribution is symmetric with respect 
to global spin reversals, we can simply fix an initial spin, 
the spin at the lower left, to have the value +1. The ori- 
entation of the remaining spins that lie along the bottom 
row of C* are then assigned sequentially first along the 
bottom row. The spins along the left column are then 
assigned. This assignment is based on the probabilities 
of domain walls separating neighboring spins in the bor- 
dering row and column. These probabilities are found by 
effectively computing the correlation functions between 
spins in the lower row and left column. Note that, in 
principle, any order of spin assignment for these outer 
border spins could be used. It is possible that numeri- 
cal stability might be improved by choosing an alternate 
order of spin assignments; we chose the nearest neighbor 
sequence for simplicity. 

Once all spins around the boundary of the sample are 
fixed, the process becomes simpler. Spin assignments are 
decided at finer scales by descending the tree recursively. 
In each subsequent step, the spins surrounding a region 
C have been fixed by prior assignment. The probabilities 
of domain wall sections crossing between the spins lying 
between two child regions A and B are computed. The 
spins between the two child regions A and B are then 
assigned by using these probabilities of relative domain 
walls. As the assignments are made, the probabilities 
for remaining parts of the separator are updated. These 
newly assigned spins then form the boundaries for the 
child regions of A and of B. These steps are shown for a 
sample spin assignment in Fig. [7] 

The iterative assignment of spins along the separators 
uses the inverse of a cluster matrix to compute correla- 
tion functions. The spins are randomly chosen accord- 
ing to these correlation functions. An essential part of 
this approach is that when a spin is fixed by such a 
choice, the inverse of the cluster matrix can be updated 
efficiently and incrementally [T2J [25] . This incremental 
update makes the sampling procedure running time for 
selecting a single spin configuration proportional to the 
time of computing the partition function (though with a 
larger prefactor) . A summary outline of the procedure is 
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Figure 7: (color online) Steps for exact sampling on a periodic Ising spin lattice of size L x X L v = 8 X 8. Empty circles indicate 
unknown spins. Arrows indicate assigned spins: f for s = +1 and J, for Sj = — 1. Separators W are drawn as dashed light 
(red) lines. The regions of the graph G* D used at each stage are shown by solid darker lines, (a) The first region C* is shown. 
The separator W (these wrapping edges are drawn as separated half edges) connects C* to itself. Given a seed spin, deciding 
which edges in W are in the matching M along the bottom row fixes the spins for step (b). (b) The spins decided in (a) fill 
the lowest row. The remaining separator along the column is to be filled in for the start of the next step, (c) In all subsequent 
steps, including this step, the boundary conditions are fixed. The region C* is separated into the lower half A* and the upper 
region B* . Choices are made for the dual edges connecting A* to B* . (d) Three edges are chosen for each of the two separators 
to fix the 3 spins for each region pair, (e) Four separators are used to set 12 spins, (f) Eight separators are used to set eight 
spins, (g) In this final stage, there are 16 separators. One edge choice is made for each, fixing the remaining undecided spins, 
(h) The final spin assignment. 
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presented in Sec. IV C| the next sections Sec. IV A and 
Sec. IV B give more details of the algorithm. 



A. Computing domain wall probabilities 

Given two regions A and B and their edge separator 
W, the domain wall probabilities are calculated using the 
inverse of Kw, the Kasteleyn matrix that incorporates 
the effects of both the values of the boundary spins sur- 
rounding A\JB and the weights on the edges in W. Note 
that boundary spins around the region A U B are taken 
to be fixed, except at the highest level. These fixed spins 
affect the weights of the intercity edges at the boundaries 
of A and B. These weights are computed for a reference 
configuration where the spins at the boundary of A and 
B are fixed to the values decided at the higher levels of 
the tree. We can continue to use all spins set to Si = +1 
for the spins interior to A and B. It follows that if a spin 
j neighbor to i in the region A is fixed to be Sj = — 1, the 
weight w(e) for the edge e crossing the bond (ij) is set to 
be exp(2/3 Jy). If the boundary spin is fixed to sj = +1, 
then the usual weight is used, w(e) = cxp(— 2/3 Jy). At 
each stage of the spin assignment, then, we recompute 
the matrices Ua and Ub using these weights that depend 
on the boundary spins for AL)B. The matrix K\y is found 
by collecting Ua and Ub into a single matrix and then 
linking the matrices using the edge weights that connect 
A and B. Matrix inversion using Pfaffian factorization is 
then used to compute the matrix . 

The inverse matrix AT^ 1 allows for the simple calcula- 
tion the probability of any given separating edge being 
part of a domain wall. These calculations Eqns. ( 16|17 ) 
use the Pfaffian analog of the Jacobi determinant iden- 
tity, which states that 



(14) 



where Uk,i is the matrix with rows and columns k and 
I removed. The notation {L/ -1 }/^ indicates the 2x2 
submatrix of U~ 1 which is built out of the intersections 
of rows and columns k and I, i.e., 



{u- l } ktl 



o u 



-u, 



k,l 



k.l 





(15) 



where k < I, so that Pf ([C/ _1 ] fc ,i) = U^\. Since the 
ratio of Pfaffians of Kasteleyn matrices Pf (if) is the ratio 
of partition functions Z, probabilities can be computed 
using the identity Eq. ( 14 ) . Let the edge on the decorated 
dual graph that separate the neighboring spins i and j 
have nodes a and b. The probability of including an edge 
e a b as part of a domain wall that separates the two spins 



is then given by the expression 

P(e ab e M) = \[K w ] ab [K 



w 



(16) 



The result Eq. ( 16 ) then follows from the probability be- 



weight of dimer configurations that don't include nodes 
k and I (i.e., Pf (Kw)k,i) divided by the total weight 
Pf (Kyy). When W is the separator for regions A and 
B with fixed boundary spins, the ab element of Ky/ is 
given by [AV] a 6 — w a b- The situation is different for the 
top level matrix C* , where W "separates" C* = A* U B* 
from itself, i.e., the edges connect boundary nodes on C* 
to each other. In this case, [A"vi/] a b = w a b + [Uc]ab- In 
addition, the probability for selecting an edge that con- 
nects C* to itself is given by a weighted sum over the 
four possible boundary dimer orientations, 



P(e ab € M) = 



E 



(r,s) 



¥f{K§ s) )[K^] 



E (r ,.)PW) 



(17) 



The equation for domain wall probabilities at the high- 
est level in the periodic lattice follows from Eq. (|5) and 
Eq. ( 14 1 and cancellations of common factors simi ar to 



those that led to the result Eq. (12 1 



ing the product of the weight of the chosen edge and the 



Once the probability of choosing an edge is computed, 
a random number y is chosen in the interval [0, 1) to de- 
cide whether to accept the addition of edge to M. If 
y < P(eij G M) the edge e is included in the sampled 
matching, otherwise it is excluded. Given the resulting 
choice, the matrix A"^ 1 is then updated by the method 
described in Sec. IIVBI We note here the contrast with 
methods for dimer covering sampling that are based on 
node separators [HUH]. In these methods, a chosen node 
that is between regions A and B was matched. The prob- 
abilities computed were the probability of choosing each 
edge that matched the chosen node, with the sum of these 
probabilities being unity. From the Jacobi identity, the 
conditional probabilities for these forced node matching 
could be found without recomputing all of the elements 
of A'^ 1 ; these probabilities are given by the Pfaffian of 
a submatrix that grows with the number of fixed nodes 
|12) . Here, instead, nodes on the boundary may or may 
not be matched, depending on whether an edge is chosen 
or not, so the same approach cannot be used. Instead, in- 
spired by the approach of Ref . [55] , we update the inverse 
matrix using the Sherman-Morrison formula. Note that 
only \W\ — 1 edges are need be chosen for each separator 
W, as the last choice of an edge is forced by consistency 
in the spin assignments (or, equivalently, parity in the 
dimer covering.) 



B. Updating K 1 using the Sherman-Morrison 
formula 

After the assignment of one spin value, the choice of 
whether the corresponding edge is included in the match- 
ing is fixed for the remainder of the calculation; all sub- 
sequent bond probabilities along the separator W must 
be computed conditioned upon this choice. This is ac- 
complished by modifying the inverse Kasteleyn matrix 
A^j} 1 for the edge separator. The Sherman-Morrison for- 
mula |25) allows for quickly recomputing the inverse of a 
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matrix when modifications of the original matrix are con- 
fined to one (or a small number) of rows and columns. 
Here, we apply this formula to set specific edges of K w to 
zero. One formulation of the Sherman- Morrison formula 
is that for any matrix A, and row vectors u and v, 



{A 



A^uv 1 A~ 
1 + v T A~ 1 i 



(18) 



In the sampling algorithm, the choice of whether to force 
the inclusion of an edge or exclude it from M modifies the 
skew-symmetric matrix Kw in two rows and columns at 
the same time. When including an edge e a b, all elements 
of Kw in rows and columns a and b are set to zero, except 
the elements [ifiy]a,& and [^iy]b,o- The matching found 
using Kw must then link a to b. In contrast, when the 
edge is excluded, these two elements [Jfyp]a,& and [i£\y]f,,a 
are set to zero, while the others in rows and columns a 
and b are kept unchanged. If a general matrix A (and 
hence its inverse A -1 ) is antisymmetric, numerical stabil- 
ity is enhanced by carrying out both row operations and 
column operations at the same time, keeping the result- 
ing matrix antisymmetric as well. Using skew-symmetry 
and applying the Sherman-Morrison formula twice gives 
the inverse of a matrix modified in two rows and columns 



as 



(A + uv T -(uv T f)~ 1 = A- 1 - 



■Fa- 1 



A^ 1 uv T A^ 1 
1 + v T A- 1 u 



(19) 



If an edge e a b is to be removed by setting K a b to zero, 
then one can set u a = 1 and Vb = —Kij, with all other 
elements of u and v being zero. Similarly, if edge ey is to 
be kept, one can also use Uk = 5 a ,k for Kronecker delta 8 
but set the vector v by Vk — —Kbkt for Vfc ^ j. 

Given that that matrix A^ 1 is computed directly only 
once, at the start of spin assignment along a separator, 
a single update procedure per Eq. [l9]may be carried out 
in 0(L 2 ) operations for a separator of L spins. This is 
faster than the 0(L 3 ) for matrix multiplication of two 
Lx L matrices because A~ 1 u and v T A~ 1 are themselves 
vectors. So performing L updates to the matrix can be 
achieved in 0(L 3 ) operations. As the nested dissection 
produces a separator with L oc yN elements, sampling 
across the separator takes 0(iV 3 / 2 ) steps. Summing this 
cost over all of the needed scales for the separators gives 
a time to sample all of the spins scaling also as C(iV 3 / 2 ). 



C. Sampling algorithm outline 

Given the connection between K w x and spin correla- 
tions and the Sherman-Morrison method for updating 
K^y as spins are chosen, the sampling algorithm for pe- 
riodic systems can now be directly described in outline 
form: 



1. Perform the same up sweep steps needed to com- 
pute the cluster matrices Ua* and Ub*- These are 
the same steps needed to compute the partition 
function Z but without the final merger. 

2. Merge Ua* and Ub* along the line of spins through 
the middle of the sample that separates the re- 
gions A* and B*. This is done by placing Ua* and 
Ub* into a larger matrix and filling in the values of 
weights w a b along this line. This gives the cluster 
matrix Uc* which is indexed by nodes along the 
bottom, top, left, and right rows of the sample. 

3. For each of the four global dimer orientations ±±, 
fill in the weights that complete the torus, using 
signs for the weights given for each choice (r, s) € 
±±. This gives four matrices Uc*,±±- Compute 
the inverse matrices Uq} using Pfaffian elimi- 
nation (Eq. j9j). 

4. For each edge e that connects the top of C* to the 
bottom of C*: 

(a) Compute the probability that the edge e is 
occupied using Eq. (17). 



(b) Apply Eq. (|19| to update U c , ± , using the 



vectors u ana v that modify Uc* so as to 
force the chosen occupation value of the cur- 
rent edge e. 

Compute two new cluster matrices, Uc* and Uc* 
which have as boundaries the left and right columns 
of C* , using the fixed values of the spins in the 
bottom row chosen in the previous step. Then use 
a modified form of Eq. |l7| ) that sums only over s, 
not both r and s, to compute probabilities for edges 
along the column at the left /right boundary of the 



sample. After selecting each edge, use Eq. ( 19 ) to 
update U^l and U^l. 

6. Use the edge choices, which give portions of rela- 
tive domain walls, to assign Ising spins around the 
border of C* . More specifically, if an edge is chosen 
to belong to the matching M, the sign of the spin 
differs on either side of the edge, while if a spin was 
not chosen, the sign of the two spins on cither side 
is the same. 

7. Sampling is now carried out recursively for the 
subregions, given these fixed boundary conditions 
around the border of the sample. This is carried 
out first for the spins lying on the central dividing 
line between A* and B* , and then for the spins ly- 
ing between their child regions, etc., until all spins 
are assigned. At each level: 

(a) Recompute the cluster matrices Ua and Ub 
for the regions A and B on either side of the 
separator. This computation uses as a refer- 
ence configuration all spins Sj = +1, except on 
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the boundary of the region A\J B, where the 
fixed boundary spins are used as the reference 
configuration. 

(b) Merge Ua and Ur using the edge weights for 
edges e £ W(A, B) between A and B to obtain 
the matrix K w , the effective Kasteleyn matrix 
for the separator W(A, B). 

(c) Compute by Pfaffian elimination. 

(d) For each of the | W\ — 1 edges e £ W, e = (i, j) 
for a node i on the boundary of A and a node 
j on the boundary of B: 

i. Compute the probability of choosing e, 
P{e) = IwijlK^ij]. 

ii. Choose whether to accept or reject the in- 
clusion of e. 

iii. Based on whether e is included or ex- 
cluded from the dimer sampling, set 
up the vectors u and v and apply 
the Sherman-Morrison formula to update 

(e) Fix the Ising spins that lie in the separator 
W using the newly computed portions of the 
domain walls. 



V. APPLICATION AND TIMING 



Jij . The sum of the Boltzmann factors at a given f3 was 
compared against the partition function Z computed for 
each sample using nested dissection and multi-precision 
arithmetic. In all cases (10 4 random samples for each 
distribution) , the partition functions were in exact agree- 
ment. We carried out tests for both bimodal and Gaus- 
sian distributions for Jy. As the support for the density 
of states is limited in the bimodal case when = ±1, 
the bimodal distribution allows for an easy exact check 
of the number of states at each energy. Setting /3 so that 
exp fi = 10 m for, say, m = 8 allows one to directly read 
off the density of states from Z written in decimal, when 
the degeneracy at all energies is less than 10 2m . 

The spin-spin correlations generated by the nested dis- 
section code were also compared with the exact enumer- 
ation results and found to be the same. In addition, to 
check our sampling code, up to 10 6 configurations were 
generated using our sampling methods for several L = 5 
samples. The temperatures used were set so that about 
90% of the configurations were in one of the ten lowest 
energy states. The temperature was set this low to have 
enough statistics to verify the Boltzmann distribution for 
the low-lying states, including the degeneracies of the bi- 
modal distribution. The distribution of energies found 
were also found to satisfy the Boltzmann distribution for 
Gaussian disorder, where there is a unique state for each 
energy. 



As the implementation of these methods into work- 
ing programs is relatively complex, we have carried out 
a number of tests of the code to confirm that it com- 
putes partition functions and correlation functions cor- 
rectly. Previous to this current version of the code, each 
of the authors has independently written a computer 
code that computes partition functions using Pfaffians. 
We confirmed that the two previous codes and the cur- 
rent partition function code [15 compute the same parti- 
tion function for samples of sizes up to size L = 256, for 
several samples at each size. We have also verified our 
code by comparison with (1) exact enumeration for small 
Ising spin glass samples of size up to 5 x 5 spins and (2) 
checking correlation functions against analytic results for 
the ferromagnetic Ising model. The following subsections 
summarize the results of tests for the pure Ising model 
and for spin glass models. Similar checks are included as 
samples in our current distribution of the partition func- 
tion code [15]. For further examples of applications of 
these particular codes, see Refs. [10| . fTT] . [26 ] [27 ] . 



A. Verification in small samples 

The checks against exact enumeration verified that the 
code produced both correct partition functions and corre- 
lation functions. A simple exact enumeration code com- 
puted the Boltzmann factor for each spin configuration S 
directly, for a given random selection of bond strengths 



B. Verification of correlations using the 
ferromagnetic model 

To check the calculation directly against an analytic 
result in larger samples, we numerically computed the 
spin-spin correlation function for the square lattice in 
ferromagnetic Ising models at the critical temperature. 
For Jij = 1 for all neighboring pairs on the square the 
lattice, the critical temperature T c satisfies T~ — j3 c — 
\ ln(\/2 + 1) for Jij = 1. The correlation function along 
the diagonals, (so,os„ in ), where the spins are now indi- 
cated by two subscripts that indicate their x and y coordi- 
nates on the lattice. The computed correlation function 
was compared with known results [551 US] ■ While this 
does not check for the effect of heterogeneities on corre- 
lation functions, it helps confirm that correlation calcu- 
lations are carried out correctly at all scales, from single 
cities up through the size of the sample. The analytic 
result |29] for an infinite sample is 



{so.qsrm) 



n 



R-l 



1 



1 

4i2 



i-R 



(20) 



Z which at large separations \i — j\ = i?\/2 3> 1 gives 

(so, sr,r) = a \i - j\- 1/4 , (21) 

with a Q — 2 1 / 12 e 3 ^ and £ is the Riemann £ function. 
The numerical results for finite-size samples are plotted 
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in Fig. [8j The rather large finite-size corrections to the 
spin-spin correlation functions are apparent, but the nu- 
merical calculation quickly converges to the exact short 
distance results of Eq. ( 20 ) and apparently converges to 



the asymptotic limit Eq. (21), giving us further confi- 



dence in the correlation function code. 
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Figure 8: Plot of correlation functions for the ferromag- 
netic Ising model at criticality. The correlation function 
(s(0, 0)s(R, R)) for diagonal spin separations (R, R) was com- 
puted using PfafHan methods for samples of size L x L spins, 
L — 8, . . . , 128, and compared with L — oo exact and asymp- 
totic results. 



C. Timings 

We conclude the review of these algorithms with a 
list of the timings and memory used, in order to com- 
pare with other implementations and algorithms. Table 
[I] gives average run times on a single core of a 2.4 GHz 
Xeon E5620 quad-core processor with 12 GB of memory. 
The GNU multiprccision arithmetic library gmp (version 
3.5.0) was utilized for high precision floating point arith- 
metic using the C++ interface gmpxx (version 4.1.0) in- 
cluded with gmp [30 . Multiprecision arithmetic was used 
for all floating point calculations, including the tempera- 
ture parameters, storing the bond strengths and weights, 
and all matrix and partial Pfaffian operations. A pair 



of custom routines were written for the logarithm and 
exponential functions. These are used only for comput- 
ing weights at the start of the computation and com- 
puting logarithms of the partition functions at the end 
of the calculation, for finding free energies. These tim- 
ings are all for periodic samples: for samples with free or 
fixed boundary conditions, approximately four times less 
memory and CPU time are needed. As long as there are 
no overflows, the running time is independent of inverse 
temperature /3, though increasing the precision increases 
the maximum value of j3 for which the calculations are 
stable. For L = 256 and bimodal disorder, floats using 
1536 bits are needed to reliably sample configurations 
for /3 = 20. Computing the partition function or using 
Gaussian disorder requires just somewhat fewer bits. 



VI. FUTURE WORK 

The goal of this paper has been to present in detail nu- 
merical methods for computing thermodynamic quanti- 
ties, computing spin-spin correlation functions, and sam- 
pling configurations for two-dimensional Ising models 
with short range (planar) interactions. The development 
and explication of these methods, which incorporates 
many ideas from previous work, emphasizes the natural 
summing over various length scales. Especially at lower 
temperatures, the near cancellations that result during 
the matrix operations require matrix operations with 
multi-precision arithmetic. The precise numerical re- 
sults obtained are a great advantage for studying thermo- 
dynamic quantities, compared with traditional Markov 
Chain Monte Carlo methods, for a broad range of prob- 
lems. We have prepared a code that should be easily 
compiled to compute partition functions for the Ising 
model with arbitrary couplings on square lattices. We 
are currently preparing implementations of the correla- 
tion function and sampling codes for distribution. The 
structure of the algorithm is also very suggestive with re- 
spect to the renormalization of couplings in random mod- 
els, which might be studied directly to look for some type 
of fixed point distribution in coarse grained couplings. 
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algorithm of Ref . |14| . This work was supported in part 
by the National Science Foundation grant DMR-1006731. 
We thank the Aspen Center for Physics, supported by 
NSF grant 1066293, for its hospitality while portions of 
this paper were written up. 
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System size 


Floating point 


CPU time, 


Peak memory, 


CPU time, 


Peak memory, 


CPU time, choose 


Peak memory, choose 


L 


precision (bits) 


compute Z 


compute Z 


correlations 


correlations 


configuration 


configuration 


16 


128 


0.21 s 




3.0 s 




0.47 s 




16 


512 


0.45 s 




3.9 s 




0.83 s 




16 


2048 


2.9 s 








4.75 s 


5.1 MB 


32 


128 


1.0 s 




18 s 


13 MB 


3.7 s 


4.9 MB 


32 
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2.2 s 




33 s 


17 MB 


6.3 s 


6.8 MB 


32 
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11 MB 
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281 MB 


348 s 


87 MB 
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146 MB 
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222 MB 


256 
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304s 


135 MB 










256 


512 


605 s 


224 MB 






3418 s 


240 MB 


256 


2048 


3946 s 


580 MB 






20268 s 


899 MB 



Table I: A listing of timings and peak memory used by the Pfaffian nested dissection algorithms for different calculations. The 
resource usage shown is for the computation for a single sample defined by bond strengths Jij. Timings and memory usage 
are displayed for computing partition functions Z, for computing correlation functions between all spins at the corners of the 
geometric dissection, and for choosing a single configuration sampled exactly from the Boltzmann distribution. The results are 
listed as a function of the size of the square sample, with N = L 2 spins, and the floating point precision used. See the text for 
a brief description of the hardware and multi-precision arithmetic libraries that were used. Blank entries in this table indicate 
where measurements were not made; dashes indicate where the matrix inverse used in the sampling method was unstable. 
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